Theory of the Raman spectrum of the rotated double-layer graphene 
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We study theoretically the Raman spectrum of the rotated double-layer graphene, consisting of 
two graphene layers rotated with respect to each other by an arbitrary angle 9. We find a relatively 
simple dependence of the Raman G peak intensity on the angle 9. On the other hand, the Raman 
2D peak position, intensity, and width show a much more complicated dependence on the angle 
9. We account for all of these effects, including dependence on the incoming photon energy, in 
good agreement with the experimental data. We find that it is sufficient to include the interaction 
between the graphene layers on the electronic degrees of freedom (resulting in the occurrence of 
van Hove singularities in the density of states). We assume that the phonon degrees of freedom 
are unaffected by the interaction between the layers. Furthermore, we decompose the Raman 2D 
peak into two components having very different linewidths; these widths are almost independent 
of the angle 9. The change in the intensity and the peak position of one of these two components 
gives insight into the dependence of the overall Raman 2D peak features as a function of the angle 
9. Furthermore, we study the influence of the coherence on the Raman signal, and we separately 
study the influence of the interaction between the layers on the electron wavefunctions and energies. 
Additionally, we show regions in the phonon spectrum giving rise to the Raman 2D peak signal. 
This work provides an insight into the interplay between the mechanical degree of freedom (angle 
9) and the electronic degrees of freedom (singularities in the density of states) in rotated double- 
layer graphene. Additionally, this work provides a way to establish experimentally the value of the 
rotation angle 9 using Raman spectroscopy measurement. This procedure becomes even more robust 
if one repeats the Raman spectroscopy measurement with a different incoming photon energy. 

PACS numbers: 78.67.Wj, 73.22. Pr, 63.22.Rc, 78.30.Na 
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I. INTRODUCTION 

The electronic band structure of a single graphene layer 
near the Fermi level consists of Dirac-cone like structure 
at the Brillouin zone edge (K point). In this work we 
study the rotated double-layer graphene which consists 
of two single-layers of graphene that are rotated with re- 
spect to each other by an arbitrary angle 9. In the special 
case when 9 = 0° the Dirac cones from the two layers are 
exactly on top of each other in reciprocal space. How- 
ever, rotation of one of the graphene layers in real space 
{9 7^ 0°) is accompanied by a corresponding rotation of 
its band structure in reciprocal space (around the ori- 
gin of the reciprocal space). Therefore, when 9^0° 
Dirac cones of the two graphene layers are no longer on 
top of each other in reciprocal space, but are separated, 
proportionally to ~ sin6'/2. Nevertheless the two Dirac 
cones are still overlapping in a small region of recipro- 
cal space in between the cones. From a perturbation 
theory argument, one would expect that the interaction 
between the Dirac cones of the two graphene layers will 
be particularly strong in region where the Dirac cones 
are overlapping. Indeed, interaction between the layers 
in the overlap region opens a hybridization gap and leads 
to van Hove singularities in the density of states of the 
rotated double-layer graphene. Since the position of the 
overlap region depends on the angle 9, we expect that 
the rotated double-layer graphene will have an interest- 
ing coupling between the mechanical degree of freedom 
(angle 9) and the electronic degrees of freedom (singu- 



larities in the density of states) . Many interesting prop- 
erties of rotated double-layer graphene arise from this 
tunability, and they have recently been attracting a lot 
of interest^"^"*^. 

In this work we study theoretically the influence of 
the angle 9 on the Raman spectrum of rotated double- 
layer graphene. Raman spectroscopy is an experimental 
technique commonly used to characterize carbon based 
materials^^. Since Raman spectroscopy uses incoming 
photons with a well defined energy, one can use this spec- 
troscopy to study selectively certain regions of the elec- 
tronic density of states. Therefore, we can expect an in- 
teresting dependence of the Raman signal of the rotated 
double-layer graphene as the angle 6 is varied. Such a de- 
pendence of the Raman signal was demonstrated in some 
recent experimental studies^'^'^^'^®"^^. 

The two most prominent Raman signals in graphene 
based systems are Raman G peak (close to 1600 cm~^) 
and Raman 2D (or G') peak (close to 2700 cm"!). The 
Raman G peak in the graphene based systems is a simpler 
process than the 2D peak since it involves creation of just 
one phonon per each scattered photon. Considering mo- 
mentum conservation and assuming a negligible momen- 
tum of the photon we conclude that the created phonon 
must occur at the Brillouin zone center. On the other 
hand, Raman 2D peak involves emission of two phonons 
per each scattered photon. In this case the momentum 
conservation implies that the two emitted phonons have 
arbitrary but opposite momenta. For this reason, the Ra- 
man 2D peak is more complex than the Raman G peak. 



In this work we study both of these Raman peaks (G and 
2D) in the rotated double-layer graphene as a function 
of angle 9. We find a relatively simple dependence of 
the Raman G peak on the angle 6. Namely, when the 
incoming photon energy is comparable to the separation 
between the van Hove singularities of the rotated double- 
layer graphene, there is a significant increase (~70) in the 
G peak intensity. On the other hand, the Raman 2D peak 
intensity, position, and width show a much more compli- 
cated dependence on the angle 9. All of these features, 
including the incoming light frequency dependence, are 
well reproduced in our calculation and agree well with 
experimental data^"^. Furthermore, these results provide 
a simple way to experimentally determine the angle of a 
rotated double-layer graphene. The angle determination 
procedure becomes even more robust if one performs Ra- 
man spectroscopy with two (or more) different incoming 
photon energies. 

We compute the Raman spectra of the rotated double- 
layer graphene using a super-cell tight-binding method. 
Additionally, for the Raman 2D peak we confirm our find- 
ings using a continuum model method. In the super-cell 
method we choose special values of 9 for which there 
exists a super-periodicity between the two graphene lay- 
ers. In the super-cell method we treat this enlarged com- 
mensurate super-cell as a unit cell of our system. On 
the other hand, in the continuum model calculation we 
rely on a simple Dirac equation description of a single- 
layer graphene and we add interaction with the other 
layer in the restricted Hilbert space. These continuum 
model calculations are less numerically demanding than 
the super-cell calculations, since they do not rely on the 
super-periodicity between the two graphene layers. How- 
ever, we expect that the super-cell tight-binding method 
is more reliable, and we find that it compares better with 
experimental data. Unless explicitly mentioned, the re- 
sults reported here refer to the super-cell tight-binding 
method. 

We provide details of both the super-cell tight-binding 
method and the continuum model method in Sec. II. In 
Sec. HI we present results of the Raman G and 2D peaks 
calculations in the rotated double-layer graphene case as 
a function of the angle 9. We also provide a detailed 
analysis of these Raman peaks in the rotated double-layer 
graphene. 



II. METHODS 

In the Raman process the incoming photon creates a 
virtual electron-hole pair which then emits (or in some 
cases, absorbs) a phonon excitation quantum (or a quan- 
tum of some other excitation). In the first-order Raman 
process, the interaction of the single incoming photon re- 
sults in the emission of a single phonon excitation. In 
the second-order Raman process two phonons are emit- 
ted for each interaction of the incoming photon. For this 
reason, measurement of the spectrum of the inelastically 



scattered outgoing photons is a sensitive probe of the 
electron and phonon degrees of freedom in the sample. 

Therefore to describe theoretically the Raman spec- 
trum of the rotated double-layer graphene, we need to 
know its electron and phonon band structures. Further- 
more, we need to evaluate the matrix element of the in- 
teraction between the electrons and light, and of the in- 
teraction between the electrons and phonons. In the re- 
mainder of this section we describe how we computed all 
of these quantities in the case of the rotated double-layer 
graphene. For simplicity, we compute the Raman inten- 
sity in the rotated double-layer graphene up to a single 
common numerical prefactor, since here we are mainly 
interested in comparing the Raman intensity of the ro- 
tated double-layer graphene to that of the single-layer 
graphene. 



A. Rotated double-layer graphene unit cell 

We start by defining the geometry of the rotated 
double-layer graphene unit cell. The single-layer 
graphene unit cell consists of two carbon atoms (A and B) 
arranged in a two-dimensional honeycomb lattice. The 
rotated double-layer graphene consists of a stack of two 
single graphene layers rotated with respect to each other 
by some angle 9. We define the angle 9 as follows. We 
start from two identical copies of single-layer graphene, 
translated along the axis perpendicular to the graphene 
plane. Next, we perform the rotation of one of the layers 
by angle 9 around the axis passing through one of the 
carbon atoms. We refer to these two layers either as the 
top and the bottom layer, or as L = 1 and L = 2 layer. 

In the case of our super-cell tight-binding method we 
work with angles 9 for which the resulting double-layer 
structure is periodic. As shown in Ref. 2 every peri- 
odic double-layer structure is characterized by a pair of 
integers p and q up to a relative translation of layers. 
Furthermore, angle 9 is related to these integers as 



cos 



3q^ - p2 
3g2 + p2 



(1) 



A few examples of the periodic double-layer structures 
are shown in Fig. 1, and they all have characteristic moire 
pattern resulting from the misalignment of two periodic 
structures. The primitive unit cell of the rotated double- 
layer graphene is indicated by black hexagons in Fig. 1. 
As can be seen from the Fig. 1, the primitive cell of the ro- 
tated double-layer graphene has a much larger area than 
the single-layer primitive cell containing only two carbon 
atoms. As shown in Ref. 2, area of the unit cell of the 
rotated double-layer graphene is N times larger than the 
single-layer unit cell area, where N is given as 



N ^ 



gcd(p, 3) 



[gcd{p + 3q,p-3q)]' 
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FIG. 1. (Color online.) A few examples of the rotated double- 
layer graphene for the special values of angle for which the 
structure is super-periodic. Carbon atoms from the fixed 
graphene layer are shown with the red dots. Blue dots in- 
dicate the carbon atoms in the layer rotated by the angle 6. 
The unit cell (super-cell) of the rotated double-layer graphene 
is indicated with black lines. Integers p and q for the each 
case are indicated on top of the each panel. Panel (a) corre- 
sponds to the double-layer graphene in which layers are not 
misaligned with respect to each other, = 0° (AA stack- 
ing). On the other hand, panel (i) corresponds to the case in 
which the angle is maximal, 9 = 60° (AB stacking). Panels 
(b) through (h) cover range of angles from 0° to 30°, specif- 
ically they are, 9.43° (b), 13.17° (c), 15.18° (d), 16.43° (e), 
17.90° (f), 21.79° (g), and 27.80° (h). 



Here gcd(a, b) is greatest common divisor of integers a 
and b. The smallest possible values of N compatible with 
Eq. 2 are 7, 13, 19, 31, 37, 43, and 49. 

An iV-fold increase of the real space primitive cell is 
accompanied with a folding of the electron and phonon 
band structure in the reciprocal space. The folded Bril- 
louin zone area is smaller by a factor of 1/N as com- 
pared to the single-layer Brillouin zone. The thick black 
hexagon in Fig. 2 shows the Brillouin zone of the rotated 
double-layer graphene for two choices of p and q, while 
the thick red and blue hexagons indicate the single-layer 
Brillouin zones of the two individual layers. As one can 
see from the figure, the rotated double-layer graphene 
Brillouin zone in these two cases needs to be repeated six 
[for case from Fig. 2(a)] or twelve [Fig. 2(b)] additional 
times in order to cover the same area as the underlying 
single-layer graphene Brillouin zone. Corresponding real- 
space super-cells for these two Brillouin zones are shown 
in panels (g) and (h) of Fig. 1. 

Throughout this paper we use a convention in which 



the wavevector from the Brillouin zone of the rotated 
double-layer graphene is denoted by k. The reciprocal 
vector of the rotated double-layer graphene is denoted 
by G. The wavevector for the Brillouin zone of the two 
single-layer graphene sheets will be denoted by k' and k", 
for the two sheets respectively. Corresponding reciprocal 
vectors of the single-layer graphene sheets we will denote 
as G' and G" . We always assume that vectors fc, k' , fe", 
G, G' , and G" are given in the same coordinate system. 



B. Electrons 

We describe the electron wavefunctions in the rotated 
double-layer graphene using a tight-binding model tak- 
ing into account interaction between the graphene lay- 
ers. Low-energy electronic excitations in graphene are 
well described by the carbon 7r-bonds. For this reason 
our model includes only one pz orbital per carbon atom 
which we will denote simply by 0(r), for the carbon atom 
at the origin. We further assume that the 0(r) orbitals at 
the two neighboring atomic sites are orthogonal to each 
other 

Using orbitals 0(r) we construct the Bloch-like tight- 
binding basis functions Xkj{f) for each fe- vector in the 
Brillouin zone and for each site j in the rotated double- 
layer graphene unit cell (super-cell) as 



;^;,^.(r)=^e^^-(«+*^)0(r-H-t,). 



(3) 



R 



Since rotated double-layer graphene consists of two 
graphene layers and the primitive unit cell of each 
graphene layer has two carbon atoms (A and B), the unit 
cell (super-cell) of the rotated double-layer graphene has 
AN carbon atoms. Therefore index j ranges from 1 to 
47V. A sum is performed over all lattice vectors iJ, while 
the coordinate of j-th orbital in the primitive unit cell is 
given by the vector tj. 

The functions Xkj{f") satisfy the periodicity require- 
ment of the Bloch theorem so we can write the ?7i-th 
electron eigenstate ip^'i'"') simply as a linear combination 
of the basis functions Xkjii"), 



i^k{r) = Y.CZ)xkAr)- 



(4) 



The band index m again ranges from 1 to 4Af, while only 
half (2A^) of these bands are assumed to be occupied for 
undoped systems (as in the single-layer graphene case). 
Next, by solving the Schrodinger equation for the elec- 
trons in the Xkj{i') basis we obtain a set of C^" coef- 
ficients at the each vector fe of interest. In order to 
construct the Schrodinger equation, we use the Slater- 
Koster parametrization from Ref. 3 fitted to the density 
functional theory calculation of the rotated double-layer 
graphene. We also rescale all tight-binding hopping pa- 
rameters from Ref. 3 by 18% to account for the GW 
computed self-energy effects^^"^^. 



1. Unfolding of the electron states 

We now describe a procedure in which one can rewrite 
(unfold) the rotated double-layer graphene electron wave- 
function in terms of the single-layer graphene basis func- 
tions. This procedure will be useful later in the computa- 
tion of the electron-phonon matrix element of the rotated 
double-layer graphene. 

In a close relation to the Eq. 3 let us now define the 
following basis functions 

a.(r) = E E e*(^+*^- V(r- R *,)• (5) 

Here, a = {L, A) is a composite index where L = 1,2 
indexes the graphene layers and X = A,B indexes the A 
and B carbon atoms. The sum over index j in Eq. 5 is 
performed over all atoms of type a, i.e. over all A or B 
carbon atoms in either first or second graphene layer. 

Functions ^fca('") respect the periodicity of the single- 
layer graphene sheet in the same way that Xkj{f) respect 
the periodicity of the rotated double-layer graphene. For 
this reason if we consider indices a from either first or 
second layer (L = 1 or 2), functions 5fca(^) become 
Bloch-like tight-binding basis functions of the first or the 
second single-layer graphene [in the same way in which 
Xfcj(r) are the basis functions of the rotated double-layer 
graphene] . 

Computing the overlap between the rotated double- 
layer graphene wavefunction ?A^('r) and the single-layer 
graphene basis function S.nai'i') gives 



p = l c? = 3 



(?««IV'r)="sc<5.,fc+GPrGc 



where we have defined quantity PJ^q^ 



as 



r>m \ ^ /^7n —iG-tj 

^kGa — 2^ ^kj^ 



(6) 



(7) 



In Eq. 6, the term 5K.^k+G equals 1 only if k and k dif- 
fer by one of the rotated double-layer reciprocal vector 
G, while Use is the total number of the super-cells in the 
entire sample. Since the |^k,q) basis is complete, Eq. 6 
implies that the electron wavefunction IV'fc') of the ro- 
tated double-layer graphene can be rewritten (unfolded) 
in terms of the basis functions \S,k+Ga) of the single-layer 
graphene (here G are the reciprocal vectors of the rotated 
double- layer graphene). Furthermore, unfolding ampli- 
tude of the rotated double-layer graphene electron wave- 
function \i]j^) in terms of the single-layer graphene basis 
function \^k+Ga) is given by quantity P^^q^ defined in 
Eq. 7. 

When preforming the unfolding procedure one needs 
to make sure that there is no double counting. For this 
reason, we always unfold states using a fixed set of G 
vectors, Q, such that the following two constraints are 
satisfied. First, vectors k + G with two different choices 
of G from the set Q never differ from each other either by 
G' or G" ( reciprocal vectors of two single-layer graphene 
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FIG. 2. (Color online.) Unfolding of the rotated double-layer 
Brillouin zone (thick black line) onto two single-layer Brillouin 
zones (red and blue). Panel (a) corresponds to p = 1, q — 3, 
6 = 21.79°, N — 7 [corresponding real-space cell is shown in 
panel (g) of Fig. 1], while panel (b) corresponds to p = 3, 
q = 7, 6 = 27.80°, A^ = 13 [corresponding real-space cell 
is shown in panel (h) of Fig. 1]. Arrows indicate G vectors 
from the set Q (see main text). Determination of the set of 
vectors Q is relatively easy for these choices of p and q, while 
it becomes more involved for some other choices since some 
elements of Q must point outside of the single-layer Brillouin 



sheets). Second, every unique single-layer wavevector k' 
or k" can be written as k + G for some G from Q and 
k from the double-layer graphene Brillouin zone. Black 
arrows on Fig. 2 indicate two examples of a set of G 
vectors G satisfying these constraints. 



C. Phonons 

In this work we assume that the interaction between 
the two layers of graphene does not affect the phonon 
band structure of the rotated double-layer graphene. 
Nevertheless, working in the Brillouin zone of the rotated 
double-layer graphene we need to take into account that 
the set of the single-layer phonons at wavevectors q + G 
are folded to the wavevector q for all G from Q. Un- 
like for the electron states, which do get affected by the 
interaction between the two graphene layers, here the un- 
folding procedure corresponds simply to the relabeling of 
states. In the folded double-layer Brillouin zone we de- 
note phonon states with the wavevector label q and the 
branch label v. On the other hand, in the unfolded single- 
layer Brillouin zone, this same phonon would be labeled 
with wavevector q + G with vector G chosen from the 
set G- Therefore labels i' and G are interchangeable for a 
unique phonon branch (corresponding for example either 
to the G or the 2D peak) of the single-layer graphene. 

In our calculations of the Raman G peak we use the 
phonon frequency of the G phonon to equal 1561 cm~^, 
as found in Ref. 26. The G phonon atomic displacement 
pattern can uniquely be determined from the representa- 
tion theory analysis of the graphene space group. 

Calculation of the Raman 2D peak is more complex 



than that of the G peak, since it involves emission of 
two phonons with arbitrary (and opposite) momentum. 
For this reason we need information about the 2D phonon 
frequencies in a relatively large region of the phonon Bril- 
louin zone close to the K-point (phonons far away from 
the K-point give negligible contribution to the Raman 
2D peak). It was shown in Ref. 27 that the Raman 2D 
peak profile relies strongly on the compensation between 
the phonon and the electron trigonal warpings, which are 
shown to be of a different sign. Therefore, the Raman 2D 
spectrum in graphene is very sensitive to the details of 
the phonon band structure. For this reason we use as 
an input to our calculations the computed 2D phonon 
frequencies from Ref. 27. Atomic displacements of 2D 
phonons are computed only at the K-points of the Bril- 
louin zone by the representation theory analysis. The 
atomic displacement pattern of the phonons near the K- 
point are approximated by the displacement pattern at 
the closest K-point. 



D. Electron-light interaction 

In the weak field approximation interaction of electrons 
with light is given by the following operator 



^iig 



light 



-A- V 



(8) 



where A is the vector potential of the electromagnetic 
field of light {h = 1). Therefore the matrix element of 
this operator between two electron states (/| and \i) up to 
a constant equals P- (/|V|i), where P is the polarization 
direction of the incoming or outgoing light. Expressing 
electron states (/| and \i) in terms of the basis functions 
(j){r) we are left with computing {(j)\'V\(j)') where {(j>\ and 
\(j)') are the tight-binding basis orbitals (/)(r) at the differ- 
ent atomic sites. We work under approximation that the 
matrix element ((/)|V|0') is exactly zero between {(J3\ and 
|(/)') that are not on the first-neighbor sites in the same 
graphene layer. Additionally, assuming pz-like symme- 
try of the (j)(r) orbitals one can easily show that under 
this approximation all matrix elements (0|V|0') can be 
determined up to a single constant prefactor. 



E. Electron-phonon interaction 

Interaction between electrons and phonons i/Pj) is de- 
scribed in terms of the deformation potential 6Vq,^. The 
deformation potential is defined as a change in the effec- 
tive potential experienced by electrons as a result of the 
phonon excitation from the z/-th phonon branch with a 
wavevector q. Similarly as in the case of an electron-light 
interaction, using symmetry and taking into account the 
interaction between the nearest neighbors, the electron- 
phonon matrix elements can be computed up to a con- 
stant prefactor for both G and 2D phonon modes. 



Since we assumed no changes to the phonon band 
structure coming from the interaction between the 
graphene layers, we compute the electron-phonon ma- 
trix element using the same electron-phonon interaction 
operator as in the case of a single-layer graphene. We 
start from the electron-phonon matrix element between 
any two rotated double-layer graphene states (fc|di and 
I''' + 9)di (dropping electron band indices). 



(fc|dii/P^|fc + q)di. 



(9) 



Next we express the rotated double-layer states in the ba- 
sis of the single-layer states using Eqs. 6 and 7 (we absorb 
coefficients P^ca ^^^'^ (• • • U and | . . .)si for simplicity). 
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J2 \f« + Q + G' 
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(10) 



Furthermore, here we also drop sum over a for simplic- 
ity. Remembering that the branch index i/ for the folded 
phonon band structure is just a relabeling of vectors G3 
from the set G we can write the electron-phonon matrix 
element as 



Y, 5](fc + Gi|.ii?P'^^.,|fe 



'q+G3 



q + G 



2ls\- 



(11) 
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Next, as we mentioned earlier, we assume that the 
electron-phonon interaction operator in the rotated 
double-layer graphene is the same as in the single-layer 
graphene. For this reason, the electron-phonon matrix 
operator conserves the crystal momentum of the single- 
layer graphene. Therefore, only one of the G3 vectors 
will give a non-zero contribution to the electron-phonon 
matrix element. 



^ Y,{k + G,UHf]^^ ^Ak- 



hG2-Gil 



q + G 



2/sl- 



(12) 
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F. Raman intensity 

Using standard perturbation technique methods^^'^^, 
one can show that the intensity of the outgoing photon 
at frequency Wout for the first-order Raman process can 
be computed as 



A(Wout) ^X! 



T^^B 



AB 



(5(a;in- Wq -Wout), (13) 



while for the second-order Raman process it is given by 



/2(l^out) ^ ^ 



qi^M 



ABC 



ABC 



(5(Win - io'Lq -l^q - ^out)- 

(14) 



Here frequency of the z^-th (/x-th) phonon branch with the 
momentum q is denoted as w^ (Wg), while the incoming 



light frequency is denoted as Wm- Furthermore, here for 
simphcity we always assume that the phonon-dependent 
terms (phonon frequencies, electron-phonon matrix ele- 
ments) appearing in the first-order Raman process are 
due to the G mode, while those appearing in the second- 
order Raman process are due to the 2D mode. Scattering 
amplitudes /^^ and I'^^(j are summed over all virtually 
excited states A, B, and C . Sum in Eqs. 13 and 14 is 
performed coherently over the electron states and inco- 
herently over the phonon states. Delta functions ensure 
the conservation of energy. 

In both Eqs. 13 and 14 we focus only on the processes 
involving emission, not absorption, of phonons, and we 
work at zero temperature. Furthermore, we are neglect- 
ing the momentum of the photon. Therefore to con- 
serve total momentum, the emitted phonon in the first- 
order Raman process must have zero momentum. In the 
second-order process momentum of one phonon (q) must 
be compensated by that of the other phonon (— q). For 
this reason, the first sum in Eq. 13 is performed over 
zero-momentum phonons, from arbitrary phonon branch 
V. Similarly, first sum in the Eq. 14 is performed over all 
pairs of phonons with momenta q and — q, from possibly 
different phonon branches v and /i. 

Scattering amplitudes I'^g and I'^^(j are most easily 
represented graphically using Feynman diagrams as in 
Figs. 3 and 4. Explicit expressions for these diagrams 
can be found in Refs. 28 and 27, here we provide as an 
example contribution from Fig. 3(a) for the first order 
Raman process, 
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FIG. 3. Feynman diagrams included in the first order Raman 
calculation for the Raman G peak. Time is increasing from 
the left to the right, photons are indicated with wavy line 
while phonon is shown with dashed line. Electrons and holes 
are drawn with arrows in the opposite direction with respect 
to time. Explicit expression for the diagram in the panel (a) 
is given in Eq. 15. 



Similarly, we also provide an explicit expression for the 
contribution of the second order Raman process from 
Fig. 4(a), 

4T„op ={k + qo\H':^^'\k + qp){kn\Hll\k + qo)- 
(fc + qp|iJP';jfcm)(fcm|iJ^is'^*|fcn)- 
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''k+q "•" '^k+q 



^in-OJ-q-e'k + q + ik-^2 



T 



(16) 



Electron bands indices in Eqs. 15 and 16 are m,n,o,p, 
while the phonon branch indices are v and /x. Electron 
eigenenergy at wavevector fc and for the band m is indi- 
cated with e^. Sum of the electron and the hole linewidth 
is given by 7, which we discuss in more detail in Sec. II G. 
Empty electron states in Eqs. 15 and 16 have the symbol 
* over their band indices. 
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FIG. 4. Feynman diagrams included in the second order Ra- 
man calculation for the Raman G peak. Conventions are as 
in Fig. 3. Explicit expression for the diagram in panel (a) is 
given in Eq. 16. 



We find that it is important to include all of the first- 
order diagrams for the Raman G peak (as shown in 
Fig. 3). For the 2D peak we include only diagrams shown 
in Fig. 4, since other permutations give much smaller 
contributions in the case of a single-layer graphene, see 
Ref. 27 for more details. 



G. Remaining parameters 

Here we discuss remaining parameters and calculation 
details used in this work. For the electron and hole 
linewidth 7 appearing in Eqs. 15 and 16 for the Ra- 
man G and 2D peak intensity we use ^ = 190 meV 
and \ = 201 meV for the 1.96 eV and 2.41 eV pho- 
ton energy calculation respectively, independent of elec- 
tron wavevector k. We choose this value to repro- 
duce the Raman G peak enhancement factor (discussed 
later, in Sec. IIIB) consistent with experiment done with 
1.96 eV incoming photon energy^'^. Nevertheless, we find 
this value to be somewhat consistent with the sum of 
linewidths coming from the electron-phonon^^ (32 and 
43 meV for 1.96 eV and 2.41 eV photon energy cal- 
culation respectively) and electron-electron interaction 
(~ 100 meV^^). Using the electron linewidth coming just 
from the electron-phonon interaction (as done in Ref. 27) 
would have resulted in a much larger Raman G peak en- 
hancement. The linewidth ^ used for the 2.41 eV in- 
coming photon energy calculation was computed from 
the value for the 1.96 eV incoming photon energy by in- 
cluding the difference in the estimated electron-phonon 
linewidths (43 meV— 32 meV=ll meV). 

To speed up the convergence of the Raman calcula- 
tion in the case of a rotated double-layer graphene we 
interpolate various electron and phonon quantities from 
a coarser reciprocal space grid onto a finer grid, and ne- 
glect contributions from the processes with a very low 
intensity. 



H. Continuum model method 

In this work, we also make use of the continuum model 
developed in Ref. 30, in order to confirm our super-cell 
tight-binding based calculation of the Raman 2D peak. 
The continuum model, as compared to the super-cell 
tight-binding method, uses an electron wavefunction in 
the restricted Hilbert space. The Brillouin zone folding 
in the super-cell tight-binding calculation implies that 
the interaction between the graphene layers introduces 
hybridization between states at A^ different wavevectors 
(state with wavevectors fc is hybridized with states fc-|-G, 
G ^ Q\ On the other hand, our continuum model as- 
sumes hybridization with only a subset of vectors G from 
Q. In particular, a state with wavelength k' in layer 
L = \ hybridizes only with the electron states in layer 
L = 2 with wavelength k' + G' . In the case of our con- 
tinuum model calculation, we consider three reciprocal 



vectors G' of layer L = 1 for which the norm \k' + G'\ 
is minimized. This approximation can be justified with 
a perturbative calculation, as in Ref. 30. 

Furthermore, as compared to the super-cell tight- 
binding calculation, in our continuum model calculation 
we are using a simple parametrization'^" of the interac- 
tion strength between the graphene layers that depends 
on only one parameter, c. Additionally, in our contin- 
uum model we neglect trigonal warping of the single-layer 
graphene band structure, and assume perfectly linear 
Dirac cone band structure parametrized with the band 
velocity vp. Following Ref. 30, we take vp = 10^ m/s 
and c = 0.11 eV. Since we neglect the trigonal warping 
effect, Raman G peak intensity vanishes in the contin- 
uum model. For this reason we use continuum model 
only to compute the Raman 2D peak. 

The electron-light matrix element, the electron-phonon 
matrix element, and the Raman intensity in the con- 
tinuum model are computed as in the super-cell based 
method. We use the same value of the electron linewidth 
in the two calculations. 



III. RESULTS AND DISCUSSION 

In this section we present results of our calculations of 
the Raman spectra in the rotated double-layer graphene. 

A. Electronic structure 

We start with a discussion of the electronic structure 
of the rotated double-layer graphene. Density of states 
for two selected angles 9 are given in Fig. 5 with thin red 
and blue lines, while that for the single-layer graphene is 
given with a thick black line. The density of states of the 
single-layer graphene in this range of energies linearly 
increases with the energy as one moves away from the 
Fermi level (Fermi level is at the zero energy in Fig. 5). 
This linear dependence of the density of states originates 
from the well known Dirac cones at the Brillouin zone 
corners of the single-layer graphene band structure. 

The two graphene layers in the rotated double-layer 
graphene are rotated with respect to each other by an 
angle 9. For this reason, the Dirac cones from each layer 
are not exactly on top of each other (in the reciprocal 
space) but are instead rotated with respect of each other 
by the angle 9. Therefore, the two Dirac cones are over- 
lapping only in a small region of the reciprocal space, and 
position of this overlap in the reciprocal space depends 
on the angle 9. In this overlap region interaction be- 
tween the two layers opens a hybridization gap, which in 
turn leads to the occurrence of prominent van Hove sin- 
gularities both in the occupied and empty states, whose 
position again depends on 9. For 9 = 6.01° case (red 
line in Fig. 5) van Hove singularities occur near ±0.5 eV, 
while for the 9 — 13.17° case (blue line in Fig. 5) van 
Hove singularities occur near ±1.0 eV. 
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FIG. 5. (Color online.) Density of states in our super-cell 
tight-binding model near the Fermi level (at the zero energy) 
for the rotated double-layer graphene at angles 6 = 6.01° 
(red) and 6 = 13.17° (blue). Two large van Hove singularities 
can be seen next to each other, with similar energy. Step-like 
singularity arises from the energy maximum or minimum (as 
a function of momentum) while the logarithmic divergence 
arises from the energy saddle point. As angle 9 is increased, 
these singularities move further away from the Fermi level 
(compare red and blue line in the figure). Thick black line 
is showing the density of states of a single-layer graphene, 
multiplied by two, so that it can be compared more easily to 
the rotated double-layer graphene case. 



FIG. 6. (Color online.) Calculated Raman G peak intensity 
as a function of angle 9 for two incoming photon energies 
[1.96 eV in black and 2.41 eV in red (gray)]. The range of 
angle 9 shown is from 0° to 30°. For range 30° < 9 < 60° we 
find almost the same Raman G peak intensity for 9 = 30° + A 
as for the 9 = 30° — A case. Intensity is measured relative to a 
single-layer graphene. We find ~ 70 fold enhancement in the 
Raman G peak intensity for 1.96 eV incoming photon energy 
near the critical angle, 10°. This enhancement shifts to the 
higher angles 9 for higher incoming photon energy (2.41 eV), 
consistent with shift in the van Hove singularity. At angles 
away from both sides of the critical angle, we find Raman G 
peak enhancement close to 2 (as would be expected in the 
limit of no interaction between the layers). 



B. Raman G peak 

As shown in the Fig. 5, the energy at which the 
van Hove singularities occur in the rotated double-layer 
graphene depends on the angle 9. For larger values of 9, 
the van Hove singularities occur further away from the 
Fermi level. In particular, for the larger value of angle 9 
the van Hove singularity of the occupied states are moved 
to lower energies, while those of the empty states are 
moved to the larger energies. When separation between 
the van Hove singularities of the empty and occupied 
states matches the incoming photon energy, we expect 
to see changes of the rotated double-layer graphene Ra- 
man spectrum. Angle 9 for which the incoming photon 
energy is close to the separation between the van Hove 
singularities we will refer to as the critical angle. 

The computed Raman G peak intensity in the rotated 
double-layer graphene is given in Fig. 6 as a function 
of angle 9 for two different incoming photon energies 
(black and red line). The Raman G peak intensity in 
Fig. 6 is given in terms of the intensity of a single-layer 
graphene. We find a ~ 70 fold enhancement of the Ra- 
man G peak intensity at angles 9 close to the critical an- 
gle, 10° (1.96 eV incoming photon energy, black line in 
Fig. 6). At angles below and above this critical angle we 
find that the Raman G peak enhancement factor is close 
to 2. Therefore, in that region of angles 6 Raman signal of 
the rotated double-layer graphene is almost the same as 
that of two independent graphene sheets. Furthermore, 



we also find that the G peak enhancement shifts to the 
higher angles 9 with higher incoming photon energy (red 
line in Fig. 6). This behavior we attribute to the shift in 
the energy of the van Hove singularity as a function of 
angle 9, as observed already in Fig. 5. 

Unlike the Raman 2D peak, the Raman G peak in 
graphene is a single phonon process and therefore its 
width and peak position depend solely on the phonon 
lifetime and frequency. We assumed in our calculation 
that the phonon lifetime and frequency are not affected 
by the interaction between the two graphene layers. For 
this reason, Raman G peak width and position are inde- 
pendent of the angle 9, in agreement with experimental 
observations in Ref. 13. 

We find a very strong dependence of the Raman G peak 
enhancement at the critical angle on the effective electron 
and hole linewidth 7 appearing in Eq. 15. Dependence 
of the Raman G peak enhancement at the critical angle 
on the value of parameter 7 is shown in Fig. 7. Dotted 
line in Fig. 7 is a fit to the function ^ 7^^. As already 
mentioned in Sec. II G due to this strong dependence of 
the Raman G peak enhancement on 7, we have chosen 
value of 7 which gives Raman G peak enhancement in 
agreement with experiment at 1.96 eV incoming photon 
energy. Nevertheless, the value of 7 we obtained is con- 
sistent with that obtained from the electron-phonon and 
electron-electron interaction estimates. 
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FIG. 7. Dependence of the Raman G peak enhancement (rel- 
ative to a single-layer) at the critical angle on the electron 
and the hole lifetime 7. The incoming photon energy in this 
calculation equals 1.96 eV. Fitted functional dependence of 
the Raman G peak enhancement (Gonh) is indicated with a 
dotted line, and equals Gcnh = 2.58 (eV^)(7/2)"^. 



1. Influence of coherence 

We find a large influence of coherence in tlie calcu- 
lation of the Raman G peak. (Similar observation was 
found in Ref. 31.) This is true both for the coherence be- 
tween different Feynman diagrams (shown in Fig. 3) and 
for the coherence between different electronic states ap- 
pearing in Eq. 13. Influence of both of these coherences 
is illustrated in Fig. 8. Figure 8 shows four different 
ways the sum given in Eq. 13 is performed. Horizontal 
axis of Fig. 8 shows the difference in electronic energies 
(AE) appearing in the energy denominator as in Eq. 15. 
The vertical axis of Fig. 8 shows value of the Raman G 
peak intensity, if the sum in Eq. 13 is performed over 
all pairs of electronic states with energy separation up to 
AE. Dotted lines in the Fig. 8 show the Raman intensity 
for the Raman G peak if the coherent sum appearing in 
Eq. 13, Eab^abI ' i^ replaced with an incoherent sum 
(X^AB I-^abI) ■ Solid lines show results for when the sum 
is performed coherently, as in Eq. 13. Additionally, dot- 
ted lines are downscaled 300 times in Fig. 8 so that they 
can be compared more easily to the coherent result. Blue 
(gray) lines in Fig. 8 represent Raman G peak intensity 
when the sum in Eq. 13 is performed only over two Feyn- 
man diagrams shown in Fig. 3 (a) and (b), while black 
lines shows results for all twelve first-order diagrams in 
Fig. 3. 



From Fig. 8 we can reach several conclusions about 
the Raman G peak in graphene. First, we find that the 
coherence in Eq. 13 between different electronic states 
leads to the suppression of the Raman G peak intensity 
by more than 300 times. Second, in order to achieve the 
fully converged result, we find that the sum in Eq. 13 has 
to be performed over electron-hole pairs separated in en- 
ergy more than the incoming photon energy Win (1.96 eV 



FIG. 8. (Color online.) Raman intensity for the G peak 
of a single-layer graphene computed in four different ways. 
Horizontal axis shows the difference in the electronic energies 
AE appearing in the energy denominators of the Feynman 
diagrams as in Eq. 15. The vertical axis shows the Raman 
G peak intensity if the sum Eq. 13 is performed only using 
electron-hole pairs separated in energy up to AE. Dotted 
lines shows the results when the sum in Eq. 13 is performed 
incoherently over electron and hole states. Both dotted lines 
are downscaled 300 times in intensity. Solid lines show the 
results when the sum is performed coherently. Blue (gray) 
lines show results when the sum is performed only over two 
Feynman diagrams shown in Fig. 3 (a) and (b), while black 
lines shows results when the sum is performed over all twelve 
Feynman diagrams. Comparing solid black line to other three 
lines, we see influence of the coherence in the electronic sum 
in Eq. 13, influence of all twelve Feynman diagrams from 
Fig. 3, and influence of performing the sum up to energies 
larger than the incoming photon energy ojin, (in this calcula- 
tion cuin =1.96 eV). 



in the case of Fig. 8) . This is especially true for the co- 
herent calculation (solid lines). Third, we find that if the 
sum in Eq. 13 is performed only up to the energies close 
to the incoming photon energy Win, that the sum is dom- 
inated by two diagrams shown in Fig. 3 (a) and (b). At 
energies larger than Wjn, other Feynman diagrams start 
to dominate [compare difference between blue (gray) and 
black lines]. 



C. Raman 2D peak 

Similarly as in the case of the Raman G peak, we ex- 
pect to see changes in the Raman 2D peak when the angle 
6 is close to the critical angle. In fact we find an even 
more complicated dependence of the 2D Raman peak on 
angle 9 than that of the Raman G peak. 

Comparing the first-order Raman calculation (as for 
the Raman G peak) given in Eq. 13 to the second-order 
Raman calculation (as for the 2D peak) given in Eq. 14 
we see that in the latter case the sum is performed over 
all phonon momenta q in the entire phonon Brillouin 
zone. Phonons at different momenta q have different fre- 
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FIG. 9. (Color online.) Super-cell tight-binding model calcu- 
lated position, intensity, and width of the Raman 2D peak. 
Black line indicates results for the incoming photon energy of 
1.96 eV while the red (gray) line shows results for the incom- 
ing photon energy of 2.41 eV. Horizontal axis gives angle d 
of the rotated double-layer graphene. The range of angle d 
shown is from 0° to 30°. For range 30° < 6 < 60° we find 
almost the same Raman G peak intensity for 6 — 30° -)- A as 
for the 9 — 30° — A case. To be consistent with Ref. 13 and 
other experimental work, fit was performed to the Lorentzian 
function. Similar results (especially for the position and the 
intensity) are obtained by a fit to the Gaussian function, see 
black line in Fig. 13. Intensity is defined as the area under 
the peak (not peak height). Width is defined as the full width 
at half of the peak maximum (FWHM). Peak intensity and 
peak position are defined relative to a single-layer graphene. 
See main text for more details. Corresponding experimental 
data (in good agreement with our calculation here) is shown 
in Ref. 13. 



quency, Wg, which in general would lead to the Gaussian- 
like spread in the Raman intensity /2('^out)j even if the 
phonon lifetime is infinite. This observation is not true 
for the Raman G peak since it involves only a single 
phonon frequency ujq, and therefore its Lorentzian- type 
width comes solely from the finite phonon lifetime, and 
its peak position is determined by ujq. 

The super-cell tight-binding method computed Raman 



FIG. 10. (Color online.) Position, intensity, and width of the 
Raman 2D peak using a more approximate method (contin- 
uum model). Conventions are the same as in Fig. 9, but the 
range of vertical scales is not the same as in Fig. 9. 



2D peak position, intensity, and width in the rotated 
double-layer graphene are given in Fig. 9. For all three 
features of the 2D peak we find a complex variation as 
a function of the angle 9, especially so near the critical 
angle (~ 10° for 1.96 eV incoming photon energy). Sim- 
ilarly as in the case of the Raman G peak, we find that 
these features shift to the larger angle 9 if the incoming 
photon energy is increased. Again, as in the case of the 
Raman G peak, this behavior is consistent with the angle 
9 dependent position of the van Hove singularities shown 
in Fig. 5. 

The position of the Raman 2D peak in Fig. 9 is in- 
dicated relative to the single-layer graphene case at the 
same incoming photon energy (since even for the single- 
layer case the Raman 2D peak position depends on 
the incoming photon energy). We find that the posi- 
tion of the Raman 2D peak of the rotated double-layer 
graphene is shifted to the larger energies with respect to 
the single-layer graphene case. The observed shift is non- 
monotonic, starting out small (r^ 4 cm"-'^) at large angles 
(> 20°). Close to the critical angle - 10° (for 1.96 eV 



11 



incoming photon energy) the shift in the peak position 
increases to ~ 14 cm"^ and is followed by a steep drop 
to ^ 4 cm^^ at about 7°. At even lower angles (< 7°) 
there is a steep rise in the 2D peak position. 

The intensity of the Raman 2D peak is somewhat less 
complicated than the peak position and the peak width. 
The Raman 2D peak intensity shows almost a step-like 
change close to the critical angle ~ 10° (for 1.96 eV in- 
coming photon energy), having an intensity comparable 
to two independent single-layers at higher angles, and 
~ 4 times smaller intensity at the smaller angles. 

The width of the Raman 2D peak at angles above 
15° (for 1.96 eV incoming photon energy) is compara- 
ble to that of a single-layer graphene, ^ 31 cm~^. At the 
smaller angles (< 15°) there is a sharp increase in the Ra- 
man 2D peak width. Additionally, close to 8° Raman 2D 
peak width suddenly jumps to 52 cm~^. Below 8° there 
is again a non-monotonic behavior of the width, starting 
with a decrease followed by a sharp increase below 6°. 

The results of the continuum model calculation of the 
Raman 2D peak in Fig. 10 show almost the same over- 
all features as the super-cell tight-binding calculations. 
The angle dependence of the peak position, intensity and 
width follow the same trends in both calculations, but 
the numerical values are somewhat different. In addi- 
tion, there are some spurious features present in Fig. 10 
that are not present in the super-cell tight-binding calcu- 
lation. For example, the Raman 2D peak intensity and 
width in the region from 6* = 5° to 6* = 15° show some 
small features not present in the super-cell tight-binding 
calculation. We expect that these differences are occur- 
ring due to the approximations introduced into the con- 
tinuum model calculation (see Sec. II H). In particular, 
the lack of the trigonal warping in the continuum model 
becomes especially important at large energy of the in- 
coming photons and for large angle 9. Additionally, the 
reduction of the Hilbert space becomes more important 
at low angles 9. 



1. Procedure to experimentally determine angle 9 

The dependence of the position, the intensity, and the 
width of the Raman 2D peak on the angle 9 provide a sim- 
ple route to experimentally determine angle 9. However, 
for some range of values of 9 the position and the width 
of the Raman 2D peak depend non-monotonically on the 
angle 9. Naively, one would expect that this would make 
it impossible to uniquely determine 9 in that range of an- 
gles. Nevertheless, combining all three properties of the 
Raman 2D peak (position, intensity, and width) make it 
easier to uniquely assign angle 9. Furthermore, combin- 
ing Raman measurements at two different incoming pho- 
ton energies gives additional way to uniquely determine 
the angle 9 even in the region where the position and the 
width of the Raman 2D peak depend non-monotonically 
on 9. For example, if one measures for the incoming pho- 
ton energy of 1.96 eV change in the Raman 2D peak po- 



sition of 8 cm"-*^, according to the black line in Fig. 9 this 
measurement can correspond to angle 6* of ~ 5°, ~ 10°, 
or ^ 15°. However if one repeats the measurement on 
the same sample with a larger incoming photon energy 
(for example 2.41 eV as shown by the red line in Fig. 9) 
and the change in the Raman 2D peak position becomes 
smaller than 8 cm~^, angle 9 can be assigned uniquely 
to ~ 10°. On the other hand, if the Raman 2D shift be- 
comes larger than 8 cm~^ then 9 is either ~ 5° or ^ 15°. 
Since these two angles are quite far apart (by construc- 
tion), other Raman 2D features like intensity or width 
can be used to determine which of the two angles should 
be assigned. Similar procedure can also be used for the 
non-monotonic dependence of the Raman 2D peak width. 



2. Decomposition into contributing phonons 



According to the Eq. 14 second order Raman process 
(as for the Raman 2D peak) can be decomposed into a 
decoherent sum of contributions coming for the pair of 
phonons (q, fi) and (— q, i') with opposite momenta q and 
possibly different phonon branches fi and v. Since the 
phonon branches of the 2D mode arise from the Brillouin 
zone folding, branch indices fi and v can be relabeled with 
the rotated double-layer graphene reciprocal vectors G 
(as discussed in Sec. II C). 

Figure 11 shows regions of the phonon Brillouin zone 
which contribute the most to the Raman 2D peak, for 
varying angle 9. Contributions from the phonon pair 
(q, ^) and (— q, v) is equally distributed among the un- 
folded vectors q+G with reciprocal vector G correspond- 
ing to both // and ly. Brillouin zone of both bottom (solid 
line) and top (dashed line) single-layer graphene are indi- 
cated with black lines. For simplicity, only contributions 
from phonons in one graphene layer are shown in Fig. 11, 
and only the region close to the Brillouin zone corner (K 
point) is shown. 

For large values of angle 9 [for example 9 — 27.80° in 
Fig. 11(h)] we find a characteristic triangular region (red) 
in the phonon Brillouin zone around the K-point with 
the largest contribution to the Raman 2D peak. Simi- 
lar behavior we find in the calculation of a single-layer 
graphene, as consistent with the decomposition found in 
Ref. 27. At angles smaller or equal to the critical angle, 
this triangular region is significantly modified. Largest 
modification we find when the K-point of the Brillouin 
zone of the top graphene layer is overlapping with the 
triangular region in the bottom layer [see for example 
Fig. 11(e)]. As shown in Fig. 12 this modification occurs 
precisely at the critical angle, at which the Dirac cones 
in the electron Brillouin zone are overlapping. 
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FIG. 11. (Color online.) Regions of the phonon Brillouin zone contributing the most to the Raman 2D peak, with characteristic 
triangular regions around the Brillouin zone K point. Region that contributes the most to the Raman 2D peak is shown in 
red color. Region of the Brillouin zone with intermediate intensity is shown in yellow, and that of zero intensity in blue. Color 
scale for each panel is scaled individually to the largest intensity for that panel, since otherwise the overall intensity of the 
panels for small angles would be too small (see Fig. 9 showing decrease of the Raman 2D peak intensity at small angles 9). 
The Brillouin zone of the bottom (solid line) and the top (dashed line) graphene layer are indicated. Contributions to the 2D 
peak of only one layer (bottom) are shown for simplicity, and we only show region of the Brillouin zone close to the K point 
(approximately the same region is indicated with dashed line in Fig. 12). The angle 9 is increasing going from the panel (a) to 
the panel (h) and it equals 4.41°(a), 7.93° (b), 8.61°(c), 9.43°(d), 10.42°(e), 11.64°(f), 13.17°(g), and 27.80° (h). Calculation is 
performed with the incoming photon energy of 1.96 eV. Large transfer of weight is seen close to the critical angle in panel (e) 
when the Brillouin zone K point of the top layer is overlapping with the triangular region in the phonon Brillouin zone. 



(a) 






FIG. 12. (Color online.) Panel (a) shows the sketch of the overlapping Dirac cones of the two graphene sheets (shown in red 
and blue). Dirac cones are centered at the Brillouin zone edge points K of the each graphene layer. Black arrows indicate the 
overlap region in which the interaction between the graphene layers introduces a hybridization gap in electron and hole states. 
Panel (b) shows the sketch of the isoenergy curves (for both layers) in the electron Brillouin zone separated in the energy by the 
amount equal to the incoming photon energy. The angle 9 is close to the critical angle. Brillouin zone of each layer is indicated 
with red and blue hexagons. Overlap region is indicated with the black arrow, as in panel (a). Panel (c) shows the sketch of 
the nesting vectors in the phonon Brillouin zone connecting two Dirac cones in the electron Brillouin zone (corresponding to 
the same graphene layer). By construction, phonon nesting vectors have opposite trigonal warping to that of the electrons and 
are twice as far away from the Brillouin zone edge point K. Approximately the same region of the phonon Brillouin zone as in 
Fig. 11 is indicated with a dashed line. 
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3. Peak substructure, two Gaussian components of the 2D 



Fit to two Gaussians 



Our calculations show that the profile of the Raman 
2D peak [/2('^out) in Eq- 14] can be well fitted with two 
Gaussians with varying position, intensity, and width of 
each Gaussian function. We find this to be true both 
for the single-layer graphene and for the rotated double- 
layer graphene. For the single-layer graphene importance 
of using two Gaussians as opposed to only one is more 
subtle. However, for the rotated double-layer graphene 
just below the critical angle, positions of these two Gaus- 
sians are somewhat apart from each other, leading to the 
more pronounced two-peak feature. Similar feature has 
been found in the experimental measurements, near the 
critical angle^^. Furthermore, these two Gaussian com- 
ponents of the Raman 2D peak behave differently as a 
function of angle 6 which will be of interest in analyz- 
ing angle 6 dependent data for the rotated double-layer 
graphene. 

First, let us analyze these two Gaussian components 
in the case of a single-layer graphene. We find that these 
Gaussian components in this case are centered around 
nearly the same frequency (difference is only 3.5 cm~^ 
at 1.96 eV incoming photon energy) and have nearly the 
same intensity. Additionally, we find that the width of 
one Gaussian component (narrow component) is 30 cni^^ 
while the width of the other Gaussian component (broad 
component) is almost two times larger, 59 cm~^. 

Figure 13 shows the position, the width, and the in- 
tensity of these two Gaussian components in the case 
of a rotated double-layer graphene (broad and narrow 
Gaussian components are shown with different color in 
Fig. 13). Data in Fig. 13 is shown for the super-cell tight- 
binding calculation, but similar results are obtained with 
the continuum model. 

Quite surprisingly, we find that the broad Gaussian 
component of the Raman 2D peak in the rotated dou- 
ble layer graphene is nearly independent of the angle 
6. There is an overall decrease in the intensity of the 
broad component below the critical angle (~ 10°) but 
the changes in the position and the width are almost 
negligible. 

For the narrow Gaussian component in the rotated 
double layer graphene we again find that its width al- 
most does not depend on the angle 9. On the other hand, 
the peak intensity and the peak position of the narrow 
component show a drastic change below the critical angle 
(~ 10°). In particular, exactly at the critical angle the 
narrow component nearly vanishes. Below the critical 
angle (5° < < 10°) the narrow component reappears 
but with significantly lower peak position (—3 cm~^ be- 
low the critical angle as compared to 18 cm~^ above the 
critical angle). At the even lower angle {0 < 5°) the 
narrow component nearly disappears once again. 

This appearance and disappearance of the narrow com- 
ponent gives an insight into the complex behavior of the 
overall position, intensity, and width of the Raman 2D 
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FIG. 13. (Color online.) Fit of the calculated Raman 2D peak 
to a single Gaussian (black line) and to two Gaussians (broad 
Gaussian component is in green and narrow in blue). Peak po- 
sition and intensity for all three lines are given relative to the 
single Gaussian fit of the 2D Raman peak in the single-layer 
graphene. Other conventions are as in the Fig. 9. Narrow 
Gaussian component for some values of angle 9 has negligible 
intensity, which makes fitting procedure ill-conditioned. For 
that range of angles, position and width of the narrow com- 
ponent are drawn with a straight dotted line. Calculation is 
performed for a single incoming photon energy, 1.96 eV. 



peak (black line in Fig. 13). For example, the overall in- 
crease in the width of the Raman 2D peak near the criti- 
cal angle (~ 10°) can be explained by the disappearance 
of the narrow Gaussian component at the same angle. 
Similarly, reappearance of the narrow component with 
lower frequency below the critical angle {5° < 6 < 10°) 
explains the overall change in the peak position of the Ra- 
man 2D peak. Additionally, reappearance of the narrow 
component at the lower frequency than the broad compo- 
nent is consistent with the experimentally observed two 
peak structure of the Raman 2D peak in the same range 
of angles 9. 

It is tempting to interpret the broad and narrow Gaus- 
sian components of the 2D peaks as coming from the 
corners of the triangular region (inner phonons, Ref. 27) 
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in Fig. 11 and from the triangular faces (outer phonons) 
respectively. Indeed, similar two-peak feature of the Ra- 
man 2D peak has been found in Ref. 27, but for signif- 
icantly larger incoming photon energies (3.8 eV). These 
two features of the Raman 2D peak were denoted as 2D"'" 
(inner) and 2D^ (outer) in Fig. 26 of Ref. 27. However, 
origin of the two peak features we find here is decidedly 
different. We demonstrate this by taking our single-layer 
graphene calculation and considering only small slices (in 
certain region of angles around the K-point) of the tri- 
angular regions in the phonon Brillouin zone either near 
the triangular corners or faces. We find in both cases 
that the two-Gaussian peak feature persists, with similar 
fitting parameters. 

Instead, we find that this two-peak structure of the 
Raman 2D peak originates from the sum over electron- 
hole pair states in Eq. 14 (not different phonon states 
as for the feature found in Ref. 27). In particular, we 
find that the electron-hole pairs which are separated by 
the energy close to the incoming light energy give rise 
to the narrow component of the 2D peak from Fig. 13, 
while the higher energy electron-hole pairs give rise to 
the broad component from Fig. 13. More specifically, for 
the incoming photon energy of 1.96 eV, we find that the 
narrow component of the 2D peak originates from the 
electron-hole pairs separated up to ^ 2.1 eV. Electron- 
hole pairs between ^ 2.1 and ~ 2.6 eV give rise to the 
broad component. 



4- Influence of interlay er interaction on electron 
wavefunctions and eigenenergies 

The tight-binding model of the rotated double-layer 
graphene used in our study is based on a Slater-Koster 
parametrization from Ref. 3. This parametrization as- 
signs a hopping term to any pair of Pz orbitals on two 
carbon atoms. These carbon atoms can either be in the 
same, or two different graphene layers. For this reason, 
we can selectively turn off interaction between the two 
graphene layers in our calculation by setting to zero all 
hopping terms between pair of carbon atoms in the dif- 
ferent graphene layers (interlayer hopping). 

The effect of allowing the electron interlayer hopping 
in our calculation is twofold. Firstly, it affects electron 
wavefunctions. The change in the electron wavefunctions 
modifies electron-light and electron-phonon matrix ele- 
ments, which in turn changes Raman intensity of both G 
and 2D peak, as given for example in the numerators of 
Eqs. 15 and 16. Secondly, interlayer hopping affects elec- 
tron eigenenergies. Electron eigenenergies in turn affect 
Raman G and 2D intensities through the denominators 
in for example Eqs. 15 and 16. 

Figure 14 shows which features of the Raman 2D peak 
can be explained solely by the influence of the interlayer 
hopping on the electron wavefunctions, and which by the 
influence on the electron eigenenergies. Dotted red (blue) 
line in Fig. 14 shows the Raman 2D peak position, inten- 
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FIG. 14. Calculated position, intensity, and width of the Ra- 
man 2D peak for the incoming photon energy of 1.96 eV. 
Dotted lines show results of the calculation in which the in- 
fluence of the electron hopping terms between two graphene 
layers affects only electron eigenenergies (red) or only electron 
wavefunctions (blue). See main text for more details. Other 
conventions are the same as in Fig. 9. 



sity, and width for the calculation in which the interlayer 
hopping is given only for the electron eigenenergies (elec- 
tron wavefunctions). Solid black line in these graphs are 
the same as in Fig. 9, showing the results of the full Ra- 
man 2D peak calculation (with interlayer hopping consid- 
ered both for electron eigenenergies and wavefunctions). 

From Fig. 14 we conclude that the position of the Ra- 
man 2D peak is almost completely determined by the 
influence of the interlayer hopping on the electron eigen- 
ergies. On the other hand, intensity of the Raman 2D 
peak is determined by the interlayer hopping influence 
on the electron wavefunctions. Finally, increase in the 
width of the Raman 2D peak at low angles is well 
described by the influence of the interlayer hopping on 
the electron wavefunctions. However, influence of the in- 
terlayer hopping on the electron wavefunctions does not 
reproduce feature in the Raman 2D peak width near the 
critical angle {^ 10°). 
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D. Limit of small and limit of large angles 

Here we discuss properties of the Raman 2D and G 
peaks of the rotated double-layer graphene in the limit 
of small (close to 0°) and large (close to 30°) angles 9. 
For the Raman G peak we find that in both limits (0° 
and 30°) intensity of the G peak is similar to that of a 
single-layer graphene (multiplied with number of layers 
in the rotated double-layer graphene, two). In fact, for 
the entire range of angles 0, except close to the critical 
angle, we find that the Raman G peak intensity is similar 
to that of a single-layer graphene (times two). 

The situation with the Raman 2D peak is again more 
complicated. Figure 15 shows calculated Raman 2D pro- 
files for the rotated double-layer graphene (black) shifted 
for clarity in the vertical direction proportionally to the 
value of the angle 9. The Raman 2D profile of the 
single-layer graphene (multiplied by two) is indicated 
with thicker red line in Fig. 15. From Fig. 15 one can see 
that the Raman 2D spectrum of the rotated double-layer 
graphene above sa 15° is already converging towards 
that of a single-layer graphene (red) . 

On the other hand, in the limit of a small angle 9 (close 
to 0°) Raman 2D peak intensity of the rotated double- 
layer graphene is significantly smaller than that at the 
larger angles, or that of the single-layer graphene. We 
find similar reduction in intensity in the case of the AB 
(blue in Fig. 15) and the AA (green in Fig. 15) stacked 
double-layer graphene. Additionally, peak position and 
width for small angles 9 are qualitatively similar to that 
of the AB and AA stacked double-layer graphene. Simi- 
larity with the AB and AA stacked double-layer graphene 
is not unexpected since the rotated double-layer graphene 
in the limit of very small angles 9 is composed of a hexag- 
onal super-periodic arrangements of AB and AA stacked 
regions. This pattern is already visible to some degree 
on Fig. 1(b) for the case oi 9 — 9.43° and is even more 
prominent at smaller angles 9. 

However, in the sharp contrast to the AB and AA 
stacked double-layer graphene, we find no prominent 
multi-peak structure in the case of the rotated double- 
layer graphene in the limit of a very small angle 9. 
Furthermore, double peak structure discussed earlier in 
Sec. HI C 3 is of a different origin, and separation in fre- 
quency between the two Gaussian components is much 
smaller. 



IV. SUMMARY AND OUTLOOK 

In this work we provided a theoretical description of 
the two most prominent Raman signals in rotated double- 
layer graphene (G peak and 2D peak). We find a rela- 
tively simple dependence of the Raman G peak intensity 
on the angle 9. On the other hand, position, intensity, 
and width of the Raman 2D peak as a function of angle 
9 is much more complex. All of our findings are in good 
agreement with available experimental data^'^. We trace 
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FIG. 15. (Color online.) Calculated Raman 2D profiles 
[/2(i^out = ijj) from Eq. 14] for the rotated double-layer 
graphene (thin black fines), the single-layer graphene mul- 
tiplied by two (thick red line), the AB stacked double-layer 
graphene (blue), and the AA stacked double-layer graphene 
(green). The rotated double-layer graphene spectra are 
shifted in the vertical direction, proportionally to the angle 
6, for clarity. Raman 2D profile of the single-layer graphene 
(red) is shifted vertically proportional to 6 = 30°. 



the origin of the complex dependence of the Raman 2D 
peak signal on the angle 9 by decomposing the Raman 
2D peak into two Gaussian components with quite dif- 
ferent widths that are nearly independent on the angle 9. 
In fact, strong dependence of the intensity and position 
of one of the components is responsible for the overall 
changes to the Raman 2D peak. 

Additionally, we discuss importance of coherence in 
the Raman G peak calculation. We analyze both coher- 
ence over the various electron-hole pairs, and coherence 
over the various Feynman diagrams contributing to the 
Raman G peak. In the case of the Raman 2D peak we 
analyze regions of the phonon Brillouin zone contribut- 
ing to the Raman signal, and explore the influence of the 
interlayer interaction on the electron wavefunctions and 
eigenenergies. 

Our study provides a way to experimentally determine 
angle 9 of the rotated double-layer graphene using only 
the Raman spectroscopy measurement. Angle determi- 
nation becomes even more robust if one repeats Raman 
spectroscopy measurement with a different incoming pho- 
ton energy, as discussed in Sec. Ill C 1. Finally, this work 
provides an insight into the coupling between the me- 
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chanical degree of freedom (angle 9) and the electronic 
degrees of freedom (singularities in the density of states) 
in the rotated double-layer graphene. We expect simi- 
lar effects to occur if even more layers of graphene are 
stacked on top of each other, or if different graphene- 
likc two-dimensional systems are stacked on top of each 
other. 



ACKNOWLEDGMENTS 

We thank Gregory Samsonidze for discussion and 
Francesco Mauri for sharing data on the calculated mono- 
layer graphene phonon band structure. This work was 
supported by the Director, Office of Science, Office of 
Basic Energy Sciences, Materials Sciences and Engineer- 
ing Division, U.S. Department of Energy under Contract 
No. DE-AC02-05CH11231 and by the National Science 
Foundation under grant No. DMRlO-1006184 which pro- 
vided for continuum model calculations. 



* sinisa@civet.berkeley.edu 

^ J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. 

Castro Neto, Phys. Rev. Lett. 99, 256802 (2007). 
^ S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. Rev. 

Lett. 101, 056803 (2008). 
^ G. Trambly de Laissardiere, D. Mayou, and L. Magaud, 

Nano Letters 10, 804 (2010). 

* E. J. Mele, Phys. Rev. B 81, 161405 (2010). 

^ R. Bistritzer and A. H. MacDonald, PNAS 108, 12233 

(2011). 
^ C. Berger, Z. Song, X. Li, X. Wu, N. Brown, C. Naud, 

D. Mayou, T. Li, J. Hass, A. N. Marchenkov, E. H. Conrad, 

P. N. First, and W. A. de Heer, Science 312, 1191 (2006). 
'' H. Schmidt, T. Liidtke, P. Barthold, E. McCann, V. I. 

Fal'ko, and R. J. Hang, Applied Physics Letters 93, 

172108 (2008). 

* J. Hicks, M. Sprinkle, K. Shepperd, F. Wang, A. Tejeda, 
A. Taleb-Ibrahimi, F. Bertran, P. Le Fevre, W. A. de Heer, 

C. Berger, and E. H. Conrad, Phys. Rev. B 83, 205403 
(2011). 

® G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Cas- 
tro Neto, A. Reina, J. Kong, and E. Y. Andrei, Nat Phys 
6, 109 (2010). 

^° A. Luican, G Li, A. Reina, J. Kong, R. R. Nair, K. S. 
Novoselov, A. K. Geim, and E. Y. Andrei, Phys. Rev. 
Lett. 106, 126802 (2011). 

" J. Hass, F. Varchon, J. E. Millan-Otoya, M. Sprinkle, 
N. Sharma, W. A. de Heer, C. Berger, P. N. First, L. Ma- 
gaud, and E. H. Conrad, Phys. Rev. Lett. 100, 125504 
(2008). 

12 Y. Wang, Z. Ni, L. Liu, Y. Liu, C. Cong, T. Yu, X. Wang, 

D. Shen, and Z. Shen, ACS Nano 4, 4074 (2010). 

13 K. Kim, S. Coh, L. Z. Tan, W. Regan, J. M. Yuk, E. Chat- 
terjee, M. F. Crommie, M. L. Cohen, S. G. Louie, and 
A. Zettl, Phys. Rev. Lett. 108, 246103 (2012). 

1^ R. W. Havener, H. Zhuang, L. Brown, R. G Hennig, and 

J. Park, Nano Letters 12, 3162 (2012). 
1^ G D. A. Jorio, R. Saito and M. S. Dresselhaus, Raman 



Spectroscopy in Graphene Related Systems (Wiley- VCH, 

Weinheim, 2011). 
^^ Z. Ni, Y. Wang, T. Yu, Y. You, and Z. Shen, Phys. Rev. 

B 77, 235403 (2008). 
1^ P. Poncharal, A. Ayari, T. Michel, and J.-L. Sauvajol, 

Phys. Rev. B 78, 113407 (2008). 
1* Z. Ni, L. Liu, Y. Wang, Z. Zheng, L.-J. Li, T. Yu, and 

Z. Shen, Phys. Rev. B 80, 125404 (2009). 
1^ A. K. Gupta, Y. Tang, V. H. Crespi, and P. C. Eklund, 

Phys. Rev. B 82, 241406 (2010). 
2° V. Carozo, C. M. Almeida, E. H. M. Ferreira, L. G. 

Canado, C. A. Achete, and A. Jorio, Nano Letters 11, 

4527 (2011). 

21 A. Righi, S. D. Costa, H. Chacham, C. Fantini, 
P. Venezuela, C. Magnuson, L. Colombo, W. S. Bacsa, 
R. S. Ruoff, and M. A. Pimenta, Phys. Rev. B 84, 241409 
(2011). 

22 S. G. Louie, Topics m Computational Materials Science, 
edited by C. Fong (World Scientific, Singapore, 1998). 

2^^ A. Griineis, C. Attaccalite, L. Wirtz, H. Shiozawa, 
R. Saito, T. Pichler, and A. Rubio, Phys. Rev. B 78, 
205425 (2008). 

2'* P. E. Trevisanutto, C. Giorgetti, L. Reining, M. Ladisa, 
and V. Olevano, Phys. Rev. Lett. 101, 226405 (2008). 

25 L. Yang, J. Deslippe, C.-H. Park, M. L. Cohen, and S. G. 
Louie, Phys. Rev. Lett. 103, 186802 (2009). 

2'' M. Mohr, J. Maultzsch, E. Dobardzic, S. Reich, 
L Milosevic, M. Damnjanovic, A. Bosak, M. Krisch, and 
C. Thomsen, Phys. Rev. B 76, 035439 (2007). 

2^^ P. Venezuela, M. Lazzeri, and F. Mauri, Phys. Rev. B 84, 
035433 (2011). 

2* M. S. Dresselhaus, A. Jorio, M. Hofmann, G. Dressel- 
haus, and R. Saito, Nano Letters 10, 751 (2010), pMID: 
20085345. 

2^ F. Jornada et al. (unpublished). 

3° R. Bistritzer and A. H. MacDonald, 108, 12233 (2011). 

31 D. M. Basko, New Journal of Physics 11, 095011 (2009). 

32 Private communication with K. Kim. 



